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The  direct  simulation  Monte  Carlo  (DSMC)  method  is  the  estabUshed  technique  for 
the  simulation  of  rarefied  gas  flows.  In  some  flows  of  engineering  interest,  such  as  occur 
for  aero-braking  spacecraft  in  the  upper  atmosphere,  DSMC  can  become  prohibitively 
expensive  in  CPU  time  because  some  regions  of  the  flow,  particularly  on  the  windward 
side  of  blunt  bodies,  become  collision  dominated.  As  an  alternative  to  using  a  hybrid 
DSMC  and  continuum  gas  solver  (Euler  or  Navier-Stokes  solver)  this  work  is  aimed 
at  making  the  particle  simulation  method  efficient  in  the  high  density  regions  of  the 
flow.  A  high  density,  infinite  collision  rate  limit  of  DSMC,  the  Equilibrium  Particle 
Simulation  method  (EPSM)  was  proposed  some  15  years  ago,  EPSM  is  developed  here 
for  the  flow  of  a  gas  consisting  of  many  different  species  of  molecules  and  is  shown  to 
be  computationally  efficient  (compared  to  DSMC)  for  high  collision  rate  flows.  It  thus 
offers  great  potential  as  part  of  a  hybrid  DSMC/EPSM  code  which  could  handle  flows 
in  the  transition  regime  between  rarefied  gas  flows  and  fully  continuum  flows.  As  a 
first  step  towards  this  goal  a  pure  EPSM  code  is  described.  The  next  step  of  combining 
DSMC  and  EPSM  is  not  attempted  here  but  should  be  straightforward.  EPSM  and 
DSMC  are  applied  to  Taylor-Couette  flow  with  Kn  =  0.02  and  0.0133  and  5^  =  3). 
Toroidal  vortices  develop  for  both  methods  but  some  differences  are  found,  as  might 
be  expected  for  the  given  flow  conditions.  EPSM  appears  to  be  less  sensitive  to  the 
sequence  of  random  numbers  used  in  the  simulation  than  is  DSMC  and  may  also  be 
more  dissipative.  The  question  of  the  origin  and  the  magnitude  of  the  dissipation  in 
EPSM  is  addressed.  It  is  suggested  that  this  analysis  is  also  relevant  to  DSMC  when 
the  usual  accuracy  requirements  on  the  cell  size  and  decoupling  time  step  are  relaxed 
in  the  interests  of  computational  eflhciency. 
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Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681-0001. 
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1  Introduction 


The  Direct  Simulation  Monte  Carlo  (DSMC)  method  [1]  [2]  is  widely  used  for  high  speed 
rarefied  gas  flows  for  which  the  continuum  equations  of  fluid  mechanics  are  expected  to  be 
invalid.  In  this  flow  regime,  the  continuum  gas  properties  such  as  density,  flow  velocity  and 
temperature  must  be  supplemented  with  information  about  the  distribution  of  molecular 
energies  amongst  the  molecules  of  the  gas.  In  the  case  that  the  only  relevant  energies  are 
those  of  translation  and  rotation  of  the  molecules,  the  distribution  function  /(v,£rot,r,t)  is 
defined  such  that,  at  a  given  point  r  in  the  gas,  at  a  given  instant  t,  the  expected  number 
density  (number /volume)  of  molecules  which  have  a  velocity  in  the  range  v  to  v  +  c?v  and 
also  have  a  value  of  rotational  energy  in  the  range  Srot  to  Crot  +  d^rot  is  given  by  nfdvdcrotr 
where  n  is  the  number  density  of  all  molecules  at  the  point. 

In  DSMC  the  molecules  of  a  gas  are  represented  by  a  very  much  smaller  number  of 
‘simulator  particles’.  The  particles  are  first  moved  in  collisionless  trajectories  for  a  small 
time  interval  At]  then  the  particle  velocities  are  subject  to  alteration  during  the  collision 
phase  of  the  calculation.  The  flow  field  is  divided  into  small  cells  and  representative  collisions 
are  calculated  amongst  the  particles  in  each  cell  while  conserving  energy  and  momentum  in 
each  collision.  The  number  of  collisions  calculated  in  each  cell  in  each  time  step  At  reflects 
the  appropriate  local  collision  rate  in  the  real  gas.  If  the  collision  rate  is  large  enough  the 
distribution  of  the  finite  sample  of  particle  velocities  in  a  cell  will  conform  (in  a  statistical 
sense)  to  the  local  Maxwell-Boltzmann  equilibrium  distribution  /o  (see  equation  (1)).  If 
the  collision  rate  is  low  there  will  be  little  re-distribution  of  molecular  velocities  and  /  will 
remain  close  to  the  value  established  in  the  cell  by  the  free  flight  part  of  the  simulation,  and 
can  thus  be  far  removed  from  Jo¬ 
in  many  blunt  body  atmospheric  re-entry  flows  the  gas  near  the  windward  surfaces  of 
a  spacecraft  can  be  at  a  high  density  and  the  continuum  fluid  dynamics  equations  should 
apply,  while  the  flow  near  the  leeward  surfaces  is  at  a  low  density  where  DSMC  should  be 
uspd  Tt  is  appealing  to  use  a  particle  method  throughout  these  simulations  rather  than  a 
hybrid  DSMC/continuum  code  in  which,  to  mention  some  apparent  difficulties,  the  inter¬ 
face  between  the  continuum  and  particle  solvers  could  be  bothersome  [3]  and  the  continuum 
solver  might  be  adversely  affected  by  the  statistical  noise  generated  by  the  DSMC  calcula¬ 
tions.  The  flow  in  the  high  density  regions  is  collision  dominated  and,  when  using  a  particle 
method,  prohibitively  expensive  amounts  of  CPU  time  can  be  spent  calculating  collisions.  A 
suggestion  by  Bird  [4]  to  limit  the  number  of  collisions  in  any  cell  in  any  time  step  to  (say) 
10  per  particle  has  sometimes  been  used.  The  justification  for  this  is  that  after  10  collisions 
per  particle  have  been  calculated  the  state  in  the  cell  will,  in  effect,  be  at  equilibrium  and 
more  collisions  cannot  improve  the  accuracy  of  the  simulation.  However,  even  this  collision 
limiting  method  can  be  prohibitively  expensive. 

The  purpose  of  this  work  is  to  describe  a  high  density  version  of  DSMC  which  will  be 
suitable  for  later  use  as  part  of  a  hybrid  particle  simulation  code.  The  high  density  direct 
simulation  method  is  called  the  Equilibrium  Particle  Simulation  Method  (EPSM)  and  was 
proposed  by  Pullin  [5]  as  an  infinite  collision  rate  limit  of  DSMC  for  a  single  species  of  particle 
only  but  it  contains  all  the  necessary  features  of  a  method  which  can  be  easily  combined 
with  DSMC.  In  a  hybrid  DSMC/EPSM  code  the  EPSM  routines  can  be  used  to  establish 
a  new  distribution  of  energies  amongst  the  particles  in  those  cells  for  which  the  theoretical 


collision  rate  would  require  more  than  (say)  10  collisions  per  particle  and  this  can  be  done 
without  calculating  any  collisions  at  all  in  those  cells.  The  two  developments  of  EPSM  which 
must  be  made  before  it  can  be  suitable  for  a  hybrid  DSMC/EPSM  code  are 

1.  rotation  energy  must  be  handled  in  a  manner  compatible  with  existing  DSMC  codes 
and 

2.  multiple  species  of  particles  must  be  allowed  for. 

These  two  issues  are  addressed  in  sections  (2)  and  (3). 

A  new  code  was  constructed  by  replacing  the  collision  routines  in  a  freely  available  DSMC 
code  with  EPSM  routines.  The  original  code,  called  DSMC2A,  was  written  by  G.  A.  Bird 
[2]  and  can  be  obtained  from  the  computer  disk  which  accompanies  ref.  [2].  DSMC2A  deals 
with  a  2D  axisymmetric  flow  in  a  rectangular  region  and  has  various  options  for  setting 
boundary  conditions  and  placing  a  surface  within  the  flow;  it  allows  for  different  species  of 
molecules,  any  of  which  may  possess  energy  of  molecular  rotation  and  vibration  as  well  as 
translational  energy.  Some  small  coding  errors  in  the  original  DSMC2A  code  were  corrected. 
These  are  listed  in  Appendix  B. 

The  new  EPSM  code  is  tested  for  Taylor- Couette  flow  and  the  results  compared  with  the 
DSMC  results  obtained  by  repeating  the  DSMC  calculations  made  at  low  Knudsen  number 
by  Stefanov  and  Cercignani  [6]  with  the  same  grid  size.  It  has  been  found  that  similar,  but 
not  identical,  vortices  will  form  using  EPSM  for  which  there  is  no  attempt  to  model  the 
collision  rate  (and  hence  the  dissipation)  correctly.  The  origin  of  the  dissipation  in  EPSM  is 
discussed  in  section  7;  that  discussion  has  some  implications  for  DSMC  when  it  is  applied 
in  the  near  continuum  regime. 


2  Vibration  energy  in  EPSM 


The  Equilibrium  Particle  Simulation  Method  (EPSM)  is  similar  to  DSMC  except  that  there 
is  no  calculation  of  molecular  collisions.  Instead,  where  DSMC  would  calculate  collisions,  it 
is  assumed  in  EPSM  that  the  collision  rate  is  always  large  enough  to  bring  about  within  a 
cell  a  distribution  of  particle  velocities  which  represents  (within  the  limitations  set  by  the 
finite  sample  size)  an  approximation  to  an  equilibrium  distribution.  For  a  single  species  of 
particle,  with  a  particle  mass  mj,  for  which  rotation  of  the  molecules  (with  two  degrees  of 
freedom)  is  active,  the  equilibrium  distribution  function  is 
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where  ki,  is  Boltzmann’s  constant.  The  distribution  of  a  single  component  of  velocity,  say 
Vx,  can  be  obtained  from  (1)  by  performing  the  integration  fodvydvzderot-  The 

result  is 
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Thus  a  single  component  of  molecular  velocity  is  normally  distributed  in  an  equilibrium  gas. 
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EPSM  bypasses  the  calculation  of  a  large  number  of  collisions  amongst  the  finite  sample 
of  particles  in  each  cell  by  selecting  the  new  velocity  components  from  a  normal  distribution 
in  such  a  way  that  the  mean  and  variance  of  the  finite  sample  conform  to  values  determined 
by  the  requirement  that  the  total  energy  and  momentum  be  conserved  in  each  cell.  This 
problem  is  different  from  the  simpler  problem  of  generating  a  finite  sample  of  velocities 
selected  from  a  parent  distribution  which  has  a  known  mean  and  variance.  That  is,  because 
of  the  finite  sample  of  simulator  particles  which  are  used  in  any  simulation  (or  any  ensemble 
of  simulations)  the  values  of  v  and  T  in  (1)  are  not  known  exactly  in  EPSM  or  DSMC; 
what  is  known  is  the  estimate  of  these  quantities  which  are  derived  from  the  finite  sample. 
Thus  the  estimate  of  the  macroscopic  fluid  velocity,  the  mass  average  particle  velocity 
is  given  by 

Vest  =  ^k'VklMtot  (3) 

fc=l 

where  N  is  the  number  of  particles  in  the  sample,  rrik  is  the  mass  of  the  k-th  particle  in  the 
sample  and 

N 

Mtot  =  Y^rnk  (4) 

fe=i 

is  the  total  mass  of  the  N  particles.  The  estimate  of  the  temperature.  Test,  is  derived  from 
the  total  energy  of  the  finite  sample  of  N  particles.  It  is  related  to  the  average  thermal 
energy  per  single  degree  of  freedom  for  a  single  particle,  ^  ^  (see  equation  (11)),  by 

^d.o.f  ~  2^bTest-  {^) 

In  the  original  application  of  EPSM  [5]  a  one  dimensional  flow  of  a  single  species  of  particle 
was  treated.  The  mathematical  derivation  of  a  procedure  for  selecting  velocity  components 
which  are  normally  distributed  and  are  constrained  to  have  a  specified  mean  and  variance 
was  developed  in  [7].  The  procedure  is  given  in  the  form  of  an  algorithm  in  [5];  the  algorithm 
is  also  given  here  in  Appendix  A  but  in  greater  detail,  particularly  for  the  special  cases  when 
the  number  of  velocities  to  be  set  is  low. 

In  the  original  application  of  EPSM  [5]  energy  of  molecular  structure,  such  as  rotation, 
was  not  assigned  to  individual  particles  in  the  simulation  but  rather  the  total  energy  of 
molecular  structure  was  calculated  by  assuming  each  degree  of  freedom  contained  exactly  the 
same  thermal  energy  as  was  contained  in  a  single  translational  degree  of  freedom^  In  current 
DSMC  codes  energy  of  molecular  rotation  is  assigned  to  each  particle,  and  for  compatibility 
with  these  codes  the  EPSM  code  developed  here  allows  for  molecular  rotation  energy  as 
described  below.  The  code  can  be  easily  generalized  to  any  number  of  fully  excited  (classical) 
degrees  of  freedom  but  the  issue  of  treating  vibration  energy  which  is  not  fully  excited  is  not 
dealt  with  here. 

It  is  convenient  to  represent  molecular  rotation  in  two  degrees  of  freedom  by  a  pseudo¬ 
velocity  having  two  components  (/ii,/«2),  where 


tin  fact  a  similar  argument  was  applied  to  the  energy  stored  in  the  two  redundant  velocity  components 
which  then  did  not  have  to  be  stored  for  the  one  dimensional  flow  case. 
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Then,  in  an  equilibrium  state,  the  distribution  function  is  given  by 


fo  (Vx,Vj„  Vj,//1,^2) 


2t:  kbT 


—mg  {v' X  +  v'y  +  v' ^ 


where  v'x  =  Vx  —  Vx,  v'y  =  Vy  —  Vy,  v'z  =  v^  —  Vz  are  the  three  components  of  thermal  velocity. 
If  an  approximation  to  the  equilibrium  distribution  (7)  can  be  established  in  a  cell  the  two 
components  of  pseudo- velocity  can  be  converted  to  a  single  value  of  rotational  energy  Srot 
which  is  stored  for  each  particle.  It  follows  from  (6)  and  (7)  that  the  distribution  of  Srot  so 
derived  will  conform  to  (1). 


3  Multiple  species 

The  EPSM  code  used  here  was  derived  by  replacing  the  collision  routines  in  DSMC2A 
(subroutines  COLLMR,  SELECT2A,  ELASTIC,  INELR  and  LBS)  with  a  subroutine  which, 
for  all  cells,  first  determines  the  mass  averaged  velocity  and  thermal  energy  of  all  the  particles 
in  the  cell.  Then  for  each  degree  of  freedom  for  each  species  of  particle  a  new  subroutine 
EQUIL  (which  executes  the  algorithm  in  Appendix  A)  is  called  in  order  to  establish  the 
new  distribution  of  energy  in  the  single  degree  of  freedom  amongst  all  the  particles  of  the 
species.  In  so  doing  the  mean  velocity  of  each  species  is  constrained  to  be  exactly  equal  to 
the  mass  averaged  velocity  and  the  random  energy  stored  in  each  degree  of  freedom  is  the 
same  for  each  species  and  for  each  degree  of  freedom.  Some  implications  of  this  are  discussed 
in  section  (4). 

The  input  to  the  EQUIL  subroutine  is  the  number  of  particle  velocities  to  be  set  Ng, 
the  required  mean  velocity  of  the  single  component  of  velocity  u,  and  e'  which  is  the  total 
specific  thermal  energy  stored  in  the  single  degree  of  freedom  divided  by  the  mass  of  one 
particle.  Thus  after  the  Ng  particle  velocities  have  been  set  (for  particles  all  having  the  same 
molecular  mass)  we  have 

1 

^  =  (8) 

5  j  —  l 

and 

=  ■  (9) 

Note  that  u  may  be  any  one  of  the  three  components  of  molecular  velocity  {vx,Vy,Vz)  or 
a  pseudo- velocity  component  representing  a  single  degree  of  freedom  for  rotation,  in  which 
case  u  =  0. 

The  required  value  of  e'  is  determined  before  each  call  to  EQUIL  from 

e' =  Nge'd.o.fjmg  (10) 

where  is  the  particle  mass  and  e'^.o./  is  the  mean  thermal  energy  per  single  degree  of 
freedom  given  by 

e'd.o.}.  =  E'  jv.  (11) 
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In  (11),  E'  is  the  total  thermal  energy  of  all  particles  in  the  cell  and  v  is  the  total  number 
of  degrees  of  freedom.  The  latter  is  given  by 

=  +  (12) 

fc=i 

where  is  the  number  of  fully  excited  (classical)  degrees  of  freedom  of  the  molecular  struc¬ 
ture;  thus  ^  =  0  for  monatomic  particles  such  as  argon  and  helium  and  ^  =  2  for  diatomic 
species  such  as  nitrogen  and  oxygen  at  room  temperatures. 

In  order  to  determine  the  total  thermal  energy  E',  the  total  translational  energy  Etr  of  the 
N  particles  in  the  cell  must  be  divided  into  the  bulk  translational  energy  and  the  random 
thermal  energy  relative  to  the  bulk  motion.  The  thermal  velocity^  of  the  k-th  particle  is 
defined  by 

v'a;  =  Vyt  -  Vest.  (13) 

It  follows  from  (3)  and  (13)  that 

E  =  E  =  E  =  0.  (14) 

fc=l  k=l  k=l 

The  total  translational  energy  of  the  N  particles  is  given  by 

Etr  =  lYl  ^  E 

^  k=l  k=l 

By  expressing  the  velocity  of  each  particle  in  terms  of  the  mass  average  velocity  and  thermal 
velocity  (15)  may  be  written  as 

■El.  =  5  (('’■!  + (“»  +  '’».)  +(”«  +  '’«)) 

=  1  E  m.  (ej  +  5“  +  vl)  .f  i  E  «  +  4  +  O 

^  k=l  ^  k=l 

N  N  N 

+  ^kv'^k  +  %  E  ^kVy,  +  Vz  E  '^kV^Zk 

k=l  k-1  kz=l 

where  use  has  been  made  of  the  (14)  and  where 

K  =  l'l2^kv'k^  (1'^) 

^  k=l 

is  the  total  translational  thermal  energy  in  the  cell.  The  total  energy  of  the  particles  is 
Etot  =  Etr  +  El^ti  that  is,  it  includes  thermal  energy  stored  in  the  rotation  of  the  particles. 
The  total  thermal  energy  {translation  +  rotation)  is  given  by 

E'  =  E;  +  E;,,  =  Etot  -  \Mtotv\  (18) 

t  Strictly  speaking,  v'  is  the  estimated  thermal  velocity  of  a  particle;  that  is,  it  is  the  velocity  relative  to 
the  simulation  estimate  of  the  true  gas  velocity  at  that  cell’s  location. 


4  Some  comments  on  EPSM  and  DSMC 


It  might  be  thought  that  another  method  of  setting  the  equilibrium  state  within  each  cell 
could  be  used.  Given  that  the  estimates  of  the  mean  velocity  and  the  temperature  in  each 
cell  are  known  from  (3)  and  (5),  a  set  of  equilibrium  particle  energies  could  be  selected 
directly  from  the  distribution  (1)  with  v  =  Vesj  and  T  =  Test,  in  exactly  the  same  way  that 
an  initial  equilibrium  state  of  a  gas  is  set  in  standard  DSMC.  However  in  that  case,  the  mean 
velocity  of  the  finite  sample  so  selected  will  not  equal  Vest  and  the  mean  thermal  energy  of 
the  finite  sample  of  velocities  will  not  correspond  to  Test]  in  other  words,  momentum  and 
energy  will  not  be  exactly  conserved  by  the  redistribution.  Nanbu’s  simulation  method  [8] 
is  another  direct  simulation  method  that  does  not  conserve  momentum  and  energy  exactly 
in  collisions,  and  this  seems  to  be  the  cause  of  the  greater  statistical  scatter  in  the  results 
compared  with  DSMC  [9]  [10].  It  seems  highly  likely  that  the  above  procedure  for  setting  an 
equilibrium  state  without  preserving  momentum  and  energy  exactly  would  be  unsatisfactory 
as  well. 

The  EPSM  procedure  does  more  than  conserve  energy  as  the  new  state  in  a  cell  is 
established.  By  setting  each  degree  of  freedom  separately  as  described  in  section  (3)  the 
final  state  is  slightly  different  from  the  state  that  would  be  established  in  DSMC  after  a 
large  number  of  collisions  had  been  calculated.  In  the  DSMC  equilibrium  state  the  mean 
velocity  of  each  species  would  be  approximately,  but  not  exactly,  the  same  as  the  mass 
averaged  velocity  and  the  different  degrees  of  freedom  would  not  contain  exactly  the  same 
thermal  energy.  There  is  some  indication  in  the  results  that  follow  that  EPSM  shows  less 
variation  than  DSMC  between  different  runs  (using  a  different  sequence  of  random  numbers) 
for  the  same  flow.  This  may  be  related  to  the  slightly  less  statistical  variation  that  is  allowed 
in  EPSM  when  the  velocity  distribution  function  is  represented  by  a  finite  sample  of  particle 
velocities. 

The  EPSM  algorithm  could  be  used  to  set  the  initial  state  in  DSMC  in  which  case  the 
initial  state  would  contain  an  exact  specified  amount  of  energy  and  each  degree  of  freedom 
would  contain  exactly  the  same  total  amount  of  thermal  energy.  It  is  possible  that  this  might 
reduce  the  statistical  scatter  in  DSMC  calculations  somewhat,  particularly  when  data  are 
obtained  by  averaging  over  an  ensemble  of  different  runs;  in  standard  DSMC  the  nominally 
identical  initial  states  contain  slightly  different  total  amounts  of  energy  and  momentum.  It 
might  also  be  possible,  and  useful,  for  something  similar  to  be  done  to  establish  incoming 
fluxes  at  flow  boundaries.  In  DSMC  a  selection  of  particles  from  an  equilibrium  gas  flows 
into  the  simulation  and  this  selection  contains  a  statistically  varying  amount  of  energy  and 
momentum  at  each  step,  even  though  the  boundary  condition  is  fixed;  it  might  be  preferable 
if  exactly  the  same  energy  and  momentum  was  carried  in  at  each  step  (while  the  individual 
particle  velocities  varied  between  steps). 

5  Computational  Efficiency  of  EPSM  vs.  DSMC 

The  expected  use  of  EPSM  as  part  of  a  hybrid  DSMC/EPSM  code  requires  that  EPSM  be 
more  efficient  than  DSMC  for  establishing  equilibrium.  In  this  section  the  computational 
efficiency  of  EPSM  is  investigated.  In  a  hybrid  code  EPSM  will  only  be  used  where  the 
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standard  DSMC  procedures  call  for  a  large  number  of  collisions  to  be  calculated  in  a  single 
time  step.  This  can  only  happen  if  the  normal  requirement  that  At  be  small  with  respect 
to  the  local  collision  time  is  relaxed,  at  least  in  some  part  of  the  flow.  If  this  criteria  is  not 
relaxed  it  follows  that  the  average  number  of  collisions  per  particle  will  be  less  than  1  and 
it  would  be  inappropriate  to  use  EPSM  in  such  a  situation  unless  it  were  known  in  advance 
that  an  equilibrium  state  existed  in  the  cell.  There  are  situations  where  large  cell  sizes  are 
used  in  regions  of  flow  where  the  flow  gradients  are  small,  and  in  that  case  the  time  step 
could  be  large  relative  to  the  local  collision  time.  This  is  particularly  likely  to  happen  for 
multi-grid  applications  of  DSMC. 

In  order  to  compare  computational  efficiency,  DSMC  and  EPSM  have  been  used  to 
simulate  a  gas  at  rest  throughout  the  region  between  two  stationary  cylinders.  The  gas 
conditions  were  the  same  as  the  initial  state  considered  in  section  (6)  for  Taylor- Couette 
flow.  The  expected  result  is  that  the  flow  remain  in  its  initial  state  (in  a  statistical  sense) 
for  all  values  of  the  time  step,  and  this  was  realized.  The  timing  results  are  shown  in  figure 
(1).  For  values  of  Aty/2RTi/ Xi  —  10,  for  which  the  number  of  collisions  per  particle  per 
time  step  is  approximately  10,  EPSM  requires  approximately  30%  of  the  CPU  time  required 
by  DSMC.  Thus  EPSM  is  an  excellent  alternative  to  DSMC  in  those  cells  in  which  many 
collisions  must  be  calculated  in  one  time  step.  It  should  also  be  noted  that  not  a  great  deal 
of  attention  has  yet  been  paid  to  considerations  of  computational  efficiency  in  EPSM;  there 
is  probably  more  room  for  improvement. 


6  EPSM  applied  to  Taylor- Couette  Flow 

In  order  to  demonstrate  the  use  of  EPSM  an  interesting  vortical  flow  (Taylor- Couette  flow) 
has  been  considered,  even  though  the  particular  conditions  considered  are  not  necessarily 
those  for  which  for  which  EPSM  will  be  cis  accurate  as  DSMC.  The  classical  case  of  Taylor- 
Couette  flow  occurs  when  a  gas  is  contained  between  two  coaxial  cylinders  of  infinite  length 
and  the  inner  cylinder  rotates  while  the  outer  cylinder  is  at  rest.  At  low  rotation  speeds,  a 
laminar  flow  is  induced  between  the  cylinders;  this  is  an  axisymmetric  version  of  the  familiar 
laminar  Couette  flow  between  two  plates  which  move  relative  to  each  other.  For  higher 
rotational  speeds  a  number  of  toroidal  vortices  will  develop.  At  even  higher  rotational  speeds 
the  vortices  can  break  down  and  the  flow  can  become  fully  turbulent  and  three  dimensional. 
Here  axisymmetric  flow  only  is  considered;  three  velocity  components  are  carried  with  the 
particles  in  the  simulation  (as  is  usual  for  DSMC)  and  the  boundary  condition  at  the  surface 
of  the  rotating  cylinder  imparts  a  tangential  velocity  to  the  particles.  Figure  (2)  shows  the 
computational  domain  and  an  example  of  a  time  averaged  EPSM  flow  with  typical  vortices. 

Stefanov  and  Cercignani  [6]  made  a  through  study  of  this  flow  for  a  moderately  rarefied 
gas  and  Bird  [2]  has  repeated  one  of  their  calculations.  Here  both  EPSM  and  DSMC  sim¬ 
ulations  have  been  performed.  For  DSMC  the  molecular  model  is  a  hard  sphere  with  the 
molecular  size  appropriate  to  Argon  at  a  temperature  of  273  K  [2].  The  ordinary  gas  con¬ 
stant  is  R  =  him  =  208.2  J  kg~^ K~^ .  Argon  has  no  rotational  or  vibrational  energy  and 
thus  the  ratio  of  specific  heats  is  7  =  5/3. 

The  inner  cylinder  radius  was  0.025  m  and  the  outer  cylinder  radius  was  0.05  m  so  that 
the  gap  between  the  cylinders  was  AR  =  0.025  m.  The  simulations  began  with  a  uniform 
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gas  at  rest  between  the  cylinders  and  the  rotation  of  the  inner  cylinder  began  impulsively. 
The  initial  gas  state  was  given  by  Ti  =  273  A',  and  two  values  of  initial  density  pi  were  used: 

1.  Pi  =  njm  =  1.716  x  10““*  kg  irT^,  Pi  =  0.130  x  10^  Pa  ,  Ai  =  5  x  10“^m  and 

2.  pi  =  riim  =  2.574  x  10”^  kg  m~^,  Pi  =  0.195  x  10®  Pa  ,  Xi  =  3.33  x  10“^m. 

The  tangential  speed  of  the  inner  cylinder  wall  was  1011.33  m.s“T  The  inner  and  outer 
cylinder  surface  temperatures  were  constant  at  =  Ti  and  diffuse  reflection  was  assumed. 
In  the  axial  direction  the  simulation  region  had  a  length  of  A  =  AAR  and  was  bounded  by 
specularly  reflecting  surfaces  at  each  end. 

For  the  given  geometry  the  flow  can  be  characterized  by  three  non-dimensional  parameters 

1.  the  Knudsen  number  based  on  the  gap  between  the  cylinders,  Kn  =  Xi/AR  =  0.02  or 
0.0133, 

2.  the  speed  ratio  =  K)/\/2i?Ti  =  3,  where  is  the  tangential  speed  of  the  inner 
cylinder  and 

3.  the  wall  to  gas  temperature  ratio  Tu,/Ti  =  1. 

Different  grids  were  used  for  the  different  Knudsen  number  flows;  in  both  cases  the  cells 
were  square  with  sides  equal  to  Xi  and  the  decoupling  interval  (time  step)  was  set  as  At  = 
lX/i/2RTi.  Variable  weighting  factors  were  not  used  so  the  number  of  simulator  particles 
remained  constant  throughout  the  simulation.  It  must  be  remembered  that  the  EPSM 
simulations  are  independent  of  the  collision  rate  (in  effect  the  collision  rate  is  infinite)  and 
hence  do  not  depend  on  the  Knudsen  number.  The  EPSM  results  may  depend  on  the  number 
of  simulator  particles  used,  the  cell  size  and  the  decoupling  interval.  Some  of  these  issues  are 
addressed  in  section  (7)  which  deals  with  the  nature  of  the  dissipation  inherent  in  EPSM. 

The  flow  can  be  visualized  with  the  aid  of  stream  traces  constructed  from  the  velocity 
field  in  the  x  -  r  plane.  Each  drawing  of  these  vortices  used  the  same  number  and  starting 
locations  of  the  stream  traces.  Each  figure  shows  the  normalized  time  i  -  tVu,/  {2Tr Ringer) 
at  the  end  of  the  sampling  period.  Except  where  explicitly  stated  in  the  caption  the  figures 
show  results  of  averaging  the  flow  over  a  single  revolution  of  the  inner  cylinder. 

Comparison  of  DSMC  and  EPSM  results  gives  an  indirect  method  of  determining  how 
closely  the  DSMC  calculations  are  approaching  the  equilibrium  limit  or  how  well  the  DSMC 
simulations  represent  the  small  departures  from  equilibrium  which  might  be  expected  for 
small  but  non-zero  Knudsen  numbers.  If  the  DSMC  results  were  no  different,  or  little 
different,  from  the  EPSM  results  it  would  mean  that  the  state  in  each  cell  was,  in  effect,  in 
local  equilibrium;  it  would  then  have  to  be  considered  whether  this  could  be  a  valid  result 
for  this  flow.  If  the  results  for  EPSM  and  DSMC  are  different  then  it  must  be  concluded 
that  EPSM  is  not  suitable  for  the  case  considered;  this  sort  of  information  will  be  useful  to 
determine  the  appropriate  regime  of  application  of  EPSM  in  any  hybrid  DSMC/EPSM  code. 
The  question  of  the  accuracy  of  the  DSMC  results  depends  on  the  usual  requirements  that 
the  cell  size  and  decoupling  time  step  be  small  relative  to  the  mean  free  path  and  collision 
time,  respectively;  these  requirements  are  approximately,  but  not  unambiguously,  satisfied 
in  these  simulations. 
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It  was  found  that  EPSM  required  about  10%  more  CPU  time  than  DSMC  in  these  simu¬ 
lations.  This  is  what  can  be  expected  according  to  figure  (1)  for  the  value  of  At  used.  This 
emphasizes  the  point  that  where  the  usual  accuracy  requirements  can  be  met  (or  approxi¬ 
mately  met),  EPSM  is,  in  theory,  not  only  inferior  (to  an  unknown  extent)  to  DSMC  but  is 
also  inefficient  (at  present)  compared  with  DSMC;  it  is  only  where  the  usual  requirements 
are  relaxed,  and  the  number  of  collisions  per  particle  per  cell  becomes  greater  than  1  that 
EPSM  is  more  efficient  than  DSMC  and  becomes  strictly  justified  on  the  grounds  that  DSMC 
approaches  the  EPSM  limit. 


6.1  Results  for  200  x  50  grid.  Kn  =  0.02 

This  grid  is  the  same  as  that  used  in  [2].  The  flow  was  sampled  at  intervals  of  3 At  so  an 
average  particle  in  the  undisturbed  gas  would  move  approximately  2  cell  widths  between 
samplings.  In  the  time  required  for  a  single  revolution  of  the  inner  cylinder  the  flow  was 
sampled  52  times. 

The  unsteady  development  of  the  flow  calculated  using  EPSM  is  shown  in  figures  (3), 
(4)  and  (5).  At  first  there  are  four  vortices  which  reduce  to  three  after  10  revolutions  and 
the  sizes  of  the  vortices  vary  with  time.  The  development  of  the  flow  using  DSMC  is  shown 
in  figures  (6),  (7)  and  (8).  DSMC  gives  a  very  different  result  from  EPSM;  there  are  five  or 
six  vortices  rather  than  three.  This  is  also  different  from  the  DSMC  result  of  Bird  [2],  which 
contained  three  vortices,  and  of  Stefanov  and  Cercignani  [6],  which  contains  four  vortices. 
A  DSMC  flow  is  not  exactly  reproducible  between  runs  and  the  same  applies  to  an  EPSM 
flow;  both  methods  depend  on  the  sequence  of  random  numbers  used  in  the  simulation. 
Apparently  DSMC  is  particularly  sensitive  to  this  effect  [11]  for  this  vortical  flow. 

These  calculations  were  repeated  for  both  EPSM  and  DSMC  using  a  different  sequence  of 
random  numbers  and  the  results  after  twenty  revolutions  are  shown  in  figures  (9)  and  (10); 
these  can  be  compared  with  those  shown  in  figures  (5)  and  (8).  Three  vortices  developed 
for  all  EPSM  calculations  made  for  this  case,  some  of  which  continued  for  sixty  revolutions. 
On  the  other  hand,  only  four  vortices  developed  for  DSMC  in  this  new  run.  The  DSMC 
calculations  used  four  sub-cells  for  each  of  the  50x200  cells  in  the  grid§  and  this  gives  an 
average  of  only  three  particles  per  sub-cell  which  seems  a  very  small  sample  from  which 
to  be  choosing  collision  partners;  this  small  number  of  particles  in  each  sub-cell  may  well 
be  the  source  of  the  larger  variation  between  runs  for  DSMC  than  for  EPSM.  The  EPSM 
procedure  was  implemented  over  the  entire  cell;  a  better  comparison  between  EPSM  and 
DSMC  could  be  made  if  EPSM  were  also  implemented  at  the  sub-cell  level,  which  has  not 
yet  been  done.  It  is  worth  noting  that  only  about  four  particles  per  cell  where  used  in 
Stefanov  and  Cercignani’s  simulations  [6]  which  also  differ  from  the  DSMC  results  here. 

6.1.1  Larger  number  of  simulator  particles 

This  case  {Kn  —  0.02,  50x200  grid)  was  re-calculated  using  twice  the  number  of  simulator 
particles  for  both  EPSM  and  DSMC.  Two  different  EPSM  results  are  shown  in  figures  (11) 
and  (12).  There  is  some  difference  between  these  results  and  each  differs  from  the  EPSM 

^Collision  partners  are  selected  from  the  same  sub-cell  while  the  collision  rate  for  the  entire  cell,  not  the 
individual  sub-cells,  is  correctly  simulated  for  the  conditions  averaged  over  the  entire  cell  [2]. 
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result  found  with  a  smaller  number  of  simulator  particles  for  which  there  were  three  rather 
than  the  five  or  six  vortices  shown  here.  For  DSMC  there  is  a  much  bigger  effect  when  the 
number  of  simulator  particles  is  doubled;  no  vortices  at  all  are  formed  (see  figures  (13)  and 

(14)). 

On  the  assumptions 


1.  that  DSMC  gives  a  better  realization  of  this  flow  than  EPSM  does  and 

2.  that  increasing  the  number  of  simulator  particles  increases  the  accuracy  of  the  simula¬ 
tion 


(two  assumptions  that  can  hardly  be  doubted)  there  is  a  trend  that  the  vortices  become 
smaller  (or  disappear)  as  the  accuracy  increases.  It  might  be  that  the  simulation  becomes  less 
dissipative  as  it  becomes  more  accurate  and  that  more  vortices  develop  for  a  less  dissipative 
flow  (in  effect  a  higher  Reynolds  number  flow)  and  eventually  the  vortices  become  unstable 
at  sufficiently  high  Reynolds  number.  However,  to  exclude  what  seems  a  remote  possibility, 
that  the  results  for  DSMC  in  this  case  correspond  to  a  laminar  flow  with  statistical  noise,  a 
detailed  statistical  analysis  along  the  lines  of  that  done  in  [6]  would  be  necessary. 

The  question  of  dissipation  in  EPSM  is  addressed  in  section  (7);  it  is  not  immediately 
clear  exactly  how  the  number  of  simulator  particles  Ntot  is  related  to  the  dissipation  in  these 
methods,  (although  for  partcile  simulation  methods  the  random  noise  ~  1  / \/Ntot  [2]). 


6.2  Results  for  300  x  75  grid.  Kn  =  0.0133 

Both  EPSM  and  DSMC  simulations  were  repeated  for  a  higher  density  flow  and  with  a 
smaller  cell  size.  The  decoupling  interval  was  such  that  Aty/2RTi/\i  remained  at  0.67.  The 
total  number  of  particles  was  twice  that  for  the  50x200  grid  and  the  average  number  of 
particles  per  cell  was  10.6,  slightly  less  than  the  minimum  on  the  50x200  grid  (11.9  particles 
per  cell).  The  DSMC  results  would  be  expected  to  change  for  this  different  physical  situation 
and  the  EPSM  results  would  also  be  expected  to  change  since,  if  the  argument  in  section 
(7)  is  correct,  the  effective  viscosity  (effective  mean  free  path)  in  EPSM  depends  on  the  cell 
size  and  decoupling  interval. 

The  unsteady  development  of  the  flow  for  EPSM  is  shown  in  figures  (15),  (16)  and  (17). 
This  can  be  compared  with  the  DSMC  flow  shown  in  figures  (18),  (19)  and  (20).  For  EPSM 
there  is  not  much  difference  at  I  =  20  between  the  vortices  in  figure  (17)  for  Kn  =  0.0133 
and  those  in  figure  (9)  for  Kn  =  0.02.  and  for  DSMC  there  is  not  much  difference  between 
figure  (20)  for  Kn  =  0.0133  and  figure  (10)  for  Kn  =  0.02.  Of  course  a  greater  difference 
can  be  found  by  comparing  with  the  low  density  runs  of  DSMC  for  which  five  vortices  (or  no 
vortices  at  all!)  formed.  If  the  purpose  of  this  work  were  to  investigate  this  Taylor-Couette 
flow  in  detail  or  to  investigate  fully  the  sensitivity  of  DSMC  in  this  application  it  would 
obviously  be  desirable  to  repeat  this  low  Knudsen  DSMC  calculation  with  a  larger  number 
of  simulator  particles;  it  could  well  be  that  the  vortices  disappear  as  they  did  for  the  higher 
Knudsen  number  case. 

While  it  is  difficult  to  draw  any  firm  conclusions  about  the  Taylor-Couette  flow  from  these 
results  it  is  interesting  to  note  that  the  DSMC  results  are  different  enough  from  the  EPSM 
results  (five  vortices  rather  than  three)  to  suggest  that  the  velocity  distribution  function  / 
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may  be  significantly  different  from  equilibrium  in  these  simulations  even  though  the  Knudsen 
number  at  0.0133  is  ‘small’.  However  it  must  be  remembered  that  the  product  S^i,Kn  is  more 
relevant  than  Kn  alone.  The  former  may  be  expressed  as 


SmKn 


A 


{2RT^) 


rj 


AR 


r  K; 


K;A 

{2RTi 


^/AR 


(19) 


and  is  thus  a  ratio  of 


1. 


a  characteristic  time  between  collisions  Uoii  =  A/  {2RTiY  and  a  characteristic  time  of 
the  macroscopic  flow  Ai?/14)  or 


2.  a  characteristic  distance  V^tcoii  traveled  (in  the  laboratory  reference  frame)  by  the 
molecules  between  collisions  and  a  characteristic  dimension  of  the  macroscopic  flow 
AR^. 


While  the  usual  definition  of  mean  free  path,  the  distance  traveled  between  collisions  mea¬ 
sured  in  a  frame  of  reference  moving  with  the  local  fluid  velocity,  makes  it  a  convenient  state 
property  that  combines  the  gas  density,  temperature  and  molecular  size,  it  is  not  necessar¬ 
ily  the  most  appropriate  microscopic  length  scale  for  high  speed  flows,  as  should  become 
apparent  after  writing  the  Boltzmann  equation  in  non-dimensional  form  (see  for  example 
[12]).  It  seems  unfortunate  that  the  term  ‘mean  free  path’  can  not  be  used  to  refer  to  the 
distance  between  collisions  measured  in  the  whatever  reference  frame  is  relevant  to  the  flow; 
if  it  could,  the  term  ‘mean  free  path’  might  be  better  applied  to 

A  =  Vcharicoll  (20) 


where 


Vchar  =  max 


[Vn  uid^  {2RTchar)'^ 


is  a  characteristic  speed  based  on  either  the  macroscopic  flow  speed  or  a  characteristic  thermal 
speed  of  the  molecules.  Equation  (20)  may  be  written  as 


A  =  max[l,  5"^]  A 


where  A  is  the  mean  free  path  as  conventionally  defined.  The  mean  collision  time  tcou  or  its 
inverse  could  serve  as  a  state  property  which  carries  information  about  the  gas  state  only, 
just  as  well  as  A  can. 


7  Effective  mean  free  path  and  dissipation  in  EPSM 

EPSM  is  an  infinite  collision  rate  limit  of  DSMC  with  a  finite  sample  of  simulator  particles. 
The  ‘infinite  particle’  limit  of  EPSM  is  known  as  the  Equilibrium  Flux  Method  (EFM)II  [5] 

^It  is  also  the  ratio  of  a  characteristic  shear  stress  in  the  flow,  niVw/AR,  and  a  characteristic  pressure, 
PiRTi  and  it  may  be  written  as  «  M^/Re,  where  M  is  the  Mach  number  and  Re  is  the  Reynolds  number. 

II  Not  an  exact  limit,  since  some  gradient  information  in  EPSM  carried  by  the  possible  non-uniform  spatial 
distribution  of  particles  in  a  cell  is  absent  from  EFM.  Pullin’s  EFM  has  has  been  rediscovered  in  [13]  and 
given  a  new  acronym  (KFVS).  Here  I  follow  the  convention  that  the  naming  rights  belong  to  the  originator 
of  the  method. 
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and  has  been  used  as  an  Euler  flow  solver  in  its  first  order  [14]  [15]  [16]  [17]  and  second 
order  [18]  [19]  [20]  formulations  and  as  the  basis  of  a  Navier-Stokes  solver  [21]  [22]  [23]. 
The  dissipation  inherent  in  EFM  has  been  demonstrated  for  a  special  case  [1-5]  and  can  be 
evaluated  more  generally  [24].  It  is  worth  reviewing  those  demonstrations  since  they  shed 
light  on  the  nature  of  the  dissipation  which  arises  in  EPSM  and,  in  some  circumstances,  in 
DSMC. 

If  the  decoupling  interval  (or  time  step  At)  in  EPSM  is  set  to  a  value  larger  than  the 
time  of  flight  across  a  typical  cell  then  a  typical  particle  will  move  to  a  new  cell  carrying 
with  it  momentum  and  energy  which  are  representative  of  the  conditions  in  the  cell  from 
which  it  came.  Next  its  momentum  and  energy  is  transmitted  to  its  new  cell  as  it  undergoes, 
in  effect,  an  infinite  number  of  collisions  with  all  the  other  particles  in  the  cell.  Clearly 
the  distance  traveled  between  collisions  in  EPSM  (and  in  DSMC)  cannot  be  less  than  the 
distance  traveled  by  a  typical  particle  in  time  At.  This  corresponds  to  a  minimum  mean  free 
path  (as  conventionally  defined  in  the  reference  frame  moving  with  the  mean  flow  velocity) 
of  approximately  {SRTjTr)^  At. 

On  the  other  hand,  consider  the  case  where  At  is  small  so  that  the  chance  that  particles 
will  cross  more  than  one  cell  in  any  one  time  step  is  negligible.  Consider  two  adjacent  cells  in 
an  EPSM  simulation  as  shown  in  figure  (21).  Let  a  local  set  of  axes  (n,p,  q)  be  attached  to 
the  interface  and  let  n  point  from  the  cell  denoted  (+)  to  the  cell  denoted  (— ).  Let  u„,Up,u, 
be  the  components  of  molecular  velocity  with  respect  to  the  local  axes.  All  particles  having 
c„  >  0  which  cross  the  interface  are  selected  from  the  equilibrium  distribution  established 
in  the  (+)  cell;  all  particles  crossing  the  interface  having  c„  <  0  are  selected  from  a  different 
equilibrium  distribution  Jq.  As  a  result  the  distribution  of  c„  of  those  particles  which  cross 
the  interface  is  discontinuous  at  c„  =  0  as  shown  in  the  figure. 

It  is  to  be  expected  that  the  fluxes  derived  from  this  (discontinuous,  non-equilibrium) 
distribution  of  velocities  of  those  particles  crossing  the  interface  will  contain  parts  corre¬ 
sponding  not  only  to  the  continuum  Euler  fluxes  (for  a  hypothetical  equilibrium  state  at  the 
interface)  but  extra  ‘dissipative’  fluxes.  Note  that  figure  (21)  shows  a  hypothetical  equilib¬ 
rium  distribution  derived  from  the  non-equilibrium  distribution;  the  full  set  of  dissipative 
fluxes  for  EPSM  (relative  to  this  hypothetical  equilibrium  distribution)  is  given  in  [24]  and 
these  results  could  be  used  to  derive  an  effective  mean  free  path  for  EPSM.  Here,  to  illustrate 
the  point,  a  simple  example  of  the  exchange  of  fluxes  between  cells  in  EPSM  is  considered. 

7.1  Fluxes  in  EPSM 

Let  Q  be  one  of  the  molecular  properties  m,  mv  or  m  +  £rot)  (mass,  momentum  and 
energy)  which  are  transported  as  the  particles  move.  The  net  flux  (/unit  area/unit  time) 
of  Q  across  the  interface  can  be  split  into  two  parts:  a  forward  flux  Fq  which  flows  in  the 
direction  of  n  and  a  backward  flux  Fq  which  flows  in  the  opposite  direction.  The  total  flux 
is 

Fq=F+  +  Fq  (21) 

where  the  two  parts  of  the  fluxes  can  be  obtained  from  standard  kinetic  theory  {e.g.  [2])  as 

/oo  roo  poo 

/  /  n'^  fo  Qvndvndvpdvq  (22) 

-OO  J  —  oo  Jo 
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and 


/oo  roo  rO 

I  I  ^  fo  Q'^nd'^ndVpdVq.  (23) 

-oo  J —oo  J —oo 

In  (22)  and  (23),  n  =  p/m  is  the  number  density  of  particles  and  m  is  the  mass  of  one 
particle. 

Now  consider  the  idealized  situation  in  figure  (22)  which  shows  adjacent  cells  in  a  section 
of  the  flow  where  there  the  gradient  dvpjdxn  0  To  evaluate  the  net  transfer  of  p- momentum 
across  the  interface  we  need  to  put  Q  =  mCp  and  perform  the  integrations  in  (22)  and  (23). 
With  the  form  of  and  fg  known  from  (1)  the  forward  and  backward  fluxes  of  p- momentum 
are  1 

FL,  =  ±  (2fir±)  "  vf  (24) 

whpfp 

W±  =  i[l±erf(s±)],  (25) 

=  ±;^exp  [- (s^)^]  .  (26) 

2X2  L  '  '  J 

and  1 

4  =  t,^/(2i{^±)^  (27) 

If  we  further  assume  there  are  no  gradients  other  than  dvpfdxn,  we  have  p+  =  p-  =  p, 

y+  —  v~  =  v'^  =  v~  =  0  and  =  T~  —  T,  and  the  net  flux  (/unit  area/unit  time)  of 

p-momentum  across  the  interface  is 

F„.,  =  p{RT/2T)i  {v*  -  v;)  .  (28) 

This  net  flux  is  the  shear  stress  r  anting  at  the  interface  and,  since  {v^  -  u")  may  be 
approximated  by  Axndvpjdxn  where  Axn  is  the  cell  size,  we  have 

T  =  p {RTI2t)'^  Axndvpidxn.  (29) 

Comparing  this  to  the  stress-strain  relationship  for  plane  Couette  flow,  r  =  p  \dv.pld\ Xn,  we 
get  the  effective  viscosity  for  EPSM  as 

p' =  {RT 12%)^  pAxn.  (30) 

Comparing  this  to  the  usual  relation  between  mean  free  path  and  viscosity 

p  ~  ^p(8i?T/x)2  X 

we  obtain  an  effective  mean  free  path  in  EPSM  of 

A'  Ax„/2.  (31) 
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Combining  the  results  for  At  >  Axjv  and  At  <  Ax/v  where  Ax  is  a  characteristic  cell 
size  and  v  is  the  magnitude  of  the  local  mean  particle  velocity  (mean  fluid  velocity)  we  can 
say  that  there  is  an  effective  mean  free  path  in  EPSM  given  by 


^EPSM 


max 


^Ax,(8i?r/7r)2  Ai 


1 

-  max 

2 


1, 


4_C 

I  S 


TV 


Ax 


(32) 


where  C  —  vAtfAx  is  the  CFL  number  and  S  =  v(  {2RT)^  is  the  local  speed  ratio.  We 
could  also  say  that  the  effective  viscosity  in  EPSM  is  given  by 


f^EPSM  ~  max 


(i?r/27r)2  pAx. 


(33) 


7.2  Implications  for  DSMC 

The  above  analysis  applies  approximately  to  DSMC  and  becomes  more  exact  as  the  distribu¬ 
tion  function  approaches  an  equilibrium  form.  Hence  the  effective  mean  free  path  in  DSMC 
would  be  given  by 

^DSMC  =  max  [0  (Ax„) ,  A] 

where  A  is  the  local  mean  free  path  for  the  molecular  collision  model  being  used.  This  is 
merely  another  way  of  stating  the  theoretical  requirement  that  the  cell  size  in  DSMC  should 
be  less  than  A.  It  is  worth  noting,  however,  that  the  effective  viscosity  in  DSMC  will  be 
given  by  an  expression  similar  to  (30)  and  not  by  the  theoretical  viscosity  of  the  molecular 
collision  model,  when  the  cell  size  is  larger  than  A,  as  may  be  the  case  when  it  is  based  on  a 
(larger)  local  characteristic  length  derived  from  flow  gradients.  In  other  words,  the  collision 
model  is  largely  irrelevant  when  the  cell  size  is  large,  and  in  that  case  EPSM  is  probably 
little  different  from  DSMC. 

It  is  obvious  that  DSMC  and  EPSM  should  give  similar  results  when  the  number  of 
collisions  per  particle  per  time  step  in  DSMC  is  large,  which  is  likely  to  be  the  case  when 
cell  sizes,  and  hence  the  decoupling  interval,  are  large.  This  does  not  necessarily  imply  an 
inadequacy  in  DSMC  or  EPSM,  since  when  flow  gradients  are  small  and  large  cell  sizes 
are  used,  the  local  flow  state  is  most  likely  to  be  very  close  to  equilibrium.  An  analogy  in 
continuum  fluid  mechanics  is  that  the  Euler  equations  (rather  than  the  full  Navier-Stokes 
equations)  can  give  reasonable  results  in  regions  of  the  flow  where  gradients  are  small. 


8  Conclusions 

The  equilibrium  particle  simulation  method  (EPSM)  [.5]  has  been  extended  to  handle  flows 
with  multiple  species  and  multiple  (classical)  degrees  of  freedom  of  molecular  structure. 
EPSM  is  now  a  suitable  method  to  be  used  as  the  high  density  part  of  a  hybrid  particle 
simulation  method  for  flows  which  contain  regions  of  rarefaction  as  well  as  regions  of  near 
continuum  flow.  EPSM  cannot  be  expected  to  be  used  as  an  alternative  to  DSMC  for  all 
flows  for  which  DSMC  might  be  used.  It  is  restricted  theoretically  to  those  cases  (which 
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may  be  only  a  small  fraction  of  the  cells  in  a  hybrid  DSMC/EPSM  simulation)  where  the 
number  of  collisions  per  particle  per  time  step  is  large  and  hence  the  distribution  function 
in  a  cell  can  be  expected  to  take  an  equilibrium  form.  In  that  case  it  has  been  found  that 
EPSM  is  very  efficient  compared  to  DSMC. 

It  was  not  the  purpose  of  this  work  to  develop  a  hybrid  DSMC/EPSM  code  but  merely 
to  demonstrate  the  use  of  EPSM.  To  this  end,  the  Taylor  Couette  flow  of  a  slightly  rarefied 
gas  {Kn  =  0.02  and  0.0133)  at  high  speed  {Syj  =  3)  has  been  calculated  using  EPSM 
and  DSMC.  Differences  between  DSMC  and  EPSM  results  were  found  and  it  follows  that 
even  at  the  low  Knudsen  number  of  0.0133  the  simulated  flow  state  is  non-equilibrium. 
What  is  more  important  than  the  Knudsen  number  alone  for  determining  the  degree  of 
departure  from  equilibrium  in  high  speed  flows  is  the  product  5'^A"n  which  is  the  ratio  of 
the  characteristic  time  between  collisions  and  a  characteristic  time  of  the  macroscopic  flow; 
in  the  flows  considered  here  this  ratio  has  a  value  of  almost  4%  which  is  apparently  enough 
to  produce  a  significant  departure  from  equilibrium  in  the  simulation.  It  doesn’t  appear  safe, 
based  on  the  variation  between  different  DSMC  results,  to  make  strong  conclusions  about 
Taylor- Couette  for  the  conditions  considered  here.  Great  computational  effort  would  need 
to  be  expended  (particularly  on  larger  numbers  of  simulator  particles  but  probably  also  on 
smaller  cell  sizes). 

EPSM  was  found  to  be  less  sensitive  to  a  change  in  the  number  of  simulator  particles  used 
than  was  DSMC.  There  was  some  indication  that  EPSM  was  more  dissipative  than  DSMC, 
possibly  because  EPSM  was  implemented  at  the  cell  level,  rather  than  the  sub-cell  level  for 
which  collisions  are  calculated  in  DSMC.  The  question  of  the  accuracy  of  EPSM  in  the  near 
continuum  regime  has  been  addressed  theoretically  by  considering  the  nature  and  magnitude 
of  the  effective  dissipation  and  mean  free  path  in  EPSM.  This  analysis  serves  to  emphasize 
some  possible  problems  in  using  particle  simulation  methods  (both  EPSM  and  DSMC)  in 
near  continuum  flows  where  the  usual  requirement  that  the  cell  size  and  decoupling  interval 
both  be  small  is  relaxed.  It  that  case  the  dissipation  in  the  simulation  methods  could  depend 
more  on  the  cell  size  than  on  the  collision  model  being  used;  in  this  respect  the  simulation 
methods  become  more  like  continuum  solution  methods  {e.g.  [15])  which  contain  an  inherent 
‘artificial’  or  ‘grid’  dissipation. 
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Appendices 


A  The  EPSM  (equilibrium)  algorithm 

The  core  of  the  EPSM  code  is  the  the  EQUIL  subroutine  which  establishes  equilibrium  values 
for  one  degree  of  freedom  for  a  finite  sample  of  particles.  The  derivation  of  the  method  is 
given  in  [7]  and  a  description  of  the  algorithm  is  given  in  [5].  A  more  detailed  pseudo-code 
description  is  given  below.  Each  call  of  TZ  returns  a  new  random  fraction  between  0  and 
1  selected  from  a  uniform  distribution.  The  input  to  the  EQUIL  subroutine  is  the  number 
of  velocity  components  to  be  set  N ^  and  the  required  values  of  u  and  the  variance  e' .  The 
output  is  the  set  of  values  uj,  j  =  1, . . .  which  are  selected  from  a  normal  distribution  and 
which  satisfy  the  constraints  u  =  Y^Uj/N  and  e'  =  ~  (see  (8)  and  (9)).  The 

algorithm  may  be  easier  to  follow  using  from  the  following  diagram  which  illustrates,  for 
particular  values  of  N,  the  order  of  generation  and  interdependence  of  the  various  terms  in 
the  algorithm. 

fV  =  10  N  =  9 

e'  u  e'  u 


Ui  I  Ui 


Interdependence  of  the  values  e,,  Vj  and  Uk  for  A’  =  10  and  N  =  9. 
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if  iV  =  2  then 

Ui  <—  u  d=  \fe'  (equal  probability  for  ±  ) 
U2  <—  2u  —  U\ 

else  if  =  3  then 


9 

•e—  2TrTZ 

V2 

4—  cos  9 

V3 

4—  V^e'sin^ 

Ui 

4—  u  —  4/202! 

U2 

*■ —  d"  (\/3^2  ~  V3)  /  V^ 

U3 

4-  U2  +  -/^Vz 

else  if  A/  =  4  then 

Ti 

62 

-e'(l-ra) 

9 

4-  2x71 

V2 

4—  •/2^cos9 

V3 

4—  4/^sm9 

Ui 

4—  u  —  v3u2/ v4u2 

U2 

4-  +  (\/4  V2  —  U3)  /  Vs 

64 

^  eTi 

V4 

4—  ±-v/2e4  (equal  probability  for  ±  ) 

U3 

4 —  U2  +  /  y/^ 

U4 

4—  U3  +  \/2  U4 

else  if . 

N  >  i  then 

if  N  is  even  then 

kmax  N  /*  maximum  index  of  Ck  */ 

for  m  =  1, . . . ,  kmaxl‘2  —  1;  step  of  1 

_ 1 _ 

Jlkmax  n-m-O.i 

end  for 

else  if  N  is  odd  then 

kmax  <—  N  —  1  /*  maximum  index  of  ek  */ 

for  m  =  1, . . . ,  kmaxl^  -  1;  step  of  1 

Tm  ^  'JZ’^’nax /2-m 

end  for 
end  if 

62  e'  —  Tf) 

9  <—  27r'R. 


V2  <—  \/2^  cos  9 
V3  <—  ^/^sm9 
ui  <—  u  —  \/N  —  1  V2ly/N 
U2  ^ +  [y/Nv2  -  IVN^ 

continued  next  page 


continue  general  case  for  >  4  .... 

P  I  /'^initialize  the  product  P  =  Ylm=i^  Tm  */ 
for  =  4, . . 

•  ?  l^Ynax  2;  step  of  2  /*last  value  of  ek,  k  =  kmax)  Is  done  later  */ 

P  ^  Tk/2-lP 

ejt  ^  e'  (l  -  Tk/2)  P 

e  ^  2%n 

Vk  +—  y/2ek  cos  6 

Vk+i  y/^  sin  6 

Uk-i  <—  Uk-2  +  (VN  +  3  —  k  Vk-i  —  y/N  -\-l  —  kv]^  jyjN  -\-2  —  k 
Uk  Uk-i  +  ^\/iV  2  —  k  Vk  —  yj N  —  k  Ufc+i)  /V N  \  —  k 

end  for 

P  ^  Tk^../2-lP 

Ck  ^  P P 

if  N  is  even  then 

VN  ^  (equal  probability  for  ±) 

else  if  N  is  odd  then 

d  e-  27rn 

'^N-l  ^V^Pk^COsd 
'^N  yJ2ek^^^  sin  9 

un-2  •*—  wat-s  +  (\/4niv-2  —  \/2vn-i  )ifm 

end  if 

un-\  ■*—  un-2  +  vm-i  —  vjv^  ly/2 

UAT  •(—  U;V-1  +  \/2  UW 

end  if  /*end  of  general  case  iV  >  4  */ 


B  Some  coding  errors  in  DSMC2A 

A  number  of  modifications  have  been  made  to  the  subroutine  M0VE2A  which  calculates 
the  movement  of  the  particles  in  the  axisymmetric  flowfield.  The  most  important  change 
concerns  the  procedures  to  be  followed  after  a  collision  with  a  reflecting  surface  has  been 
detected.  The  particle  is  moved  to  the  collision  point  by  calling  subroutine  AIFR.  The  input 
to  AIFR  is  the  initial  radial  coordinate  r  (which  corresponds  to  the  initial  y  coordinate  in 
3D  space),  the  initial  radial  and  tangential  velocity  components  (which  correspond  to  the  Vy 
and  velocity  components  in  3D  space),  and  the  position  coordinate  changes  Ay  and  A2: 
which  occur  during  the  time  of  flight  up  until  collision  with  the  surface.  The  output  from 
AIFR  is  the  new  radial  coordinate  and  the  new  radial  and  tangential  velocity  components 
(there  is  a  continuous  alteration  of  the  radial  and  tangential  coordinate  directions  as  the 
particle  moves  trough  3D  space  with  constant  velocity). 

The  first  apparent  error  in  M0VE2A  occurs  in  the  call  to  AIFR  when  a  collision  is 
detected  with  a  surface  which  is  normal  to  the  x  axis;  AIFR  is  called  with  increments  Ay 
and  Az  which  are  appropriate  to  the  time  of  flight  after  rather  than  before  the  collision.  The 
second,  and  more  serious  error,  concerns  the  checking  of  the  weighting  factors  after  surface 
collisions  to  determine  whether  the  particle  (numbered  N)  should  be  deleted  after  collision 
or  duplicate  particles  should  be  created.  In  that  the  particle  numbered  N  is  moved  for  the 
time  of  flight  remaining  after  collision  on  return  from  WEIGHT,  the  code  implicitly  assumes 
that  the  particle  was  not  deleted.  If  in  fact  the  particle  was  deleted  inside  WEIGHT  the 
value  of  N  will  have  been  decremented  by  1  and  so  a  lower  numbered  particle,  which  will 
have  been  moved  previously  in  the  given  time  step,  will  be  moved  again. 

This,  perhaps  small,  error  can  lead  to  a  more  serious  effect.  If  the  first  numbered  par¬ 
ticle  A^  =  1  happens  to  collide  with  a  reflecting  surface,  and  happens  to  be  deleted  inside 
WEIGHT,  N  will  have  the  value  0  after  the  return  from  WEIGHT.  The  code  will  attempt 
to  move  the  non-existent  particle  N  =  G]  this  leads  to  references  to  memory  locations  which 
are  ‘out  of  bounds’  of  various  arrays.  To  avoid  this  some  extra  few  lines  of  code  have  been 
added  which  detect  whether  the  particle  was  deleted  inside  WEIGHT;  if  so  the  code  moves 
on  to  the  next  particle** 

Another  very  small  error  in  the  conservation  of  mass  was  detected.  It  was  noticed  that 
for  the  closed  flow  volume  between  the  two  rotating  cylinders  of  the  Taylor- Couette  flow 
the  total  number  of  particles  in  the  simulation  was  not  constant,  even  though  no  weighting 
factors  were  used;  there  was  a  tiny  percentage  of  simulator  particles  w'hich  were  removed 
from  the  simulation.  To  overcome  this,  a  number  of  small  changes  were  made  to  the  surface 
collision  detection  logic.  These  changes  were  designed  to  catch  some  singular  cases  when 
the  particle  happened  to  pass  exactly  (to  machine  precision)  through  the  end  point  of  a 
reflecting  surface  and  when  the  particle  finished  its  flight  exactly  on  the  surface.  A  similar 
problem  occured  when  a  particle  would  finish  its  flight  just  beyond  a  surface;  it  one  part 
of  the  code  this  would  be  detected  as  a  valid  collision  (which  it  is)  but  in  a  subsequent 
part  of  the  code  (because  of  rounding  errors)  this  case  was  rejected  as  a  collision  and  the 

**The  general  axisymmetric  code  in  use  at  Imperial  College  for  many  years  [25]  [26],  which  was  derived 
from  an  early  version  of  the  codes  now  designated  G2  in  [2],  checked  weighting  factors  only  when  the  final 
position  of  the  particle  was  found.  Implementing  this  strategy  in  the  DSMC2A  code  is  another  possible  way 
of  avoiding  the  problem. 
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particle  was  left  in  its  position  just  beyond  the  surface.  Since  the  surface  lay  along  the 
boundary  of  the  simulation  region  the  particle  was  subsequently  deleted.  A  change  to  8- 
byte  (double  precision)  arithmetic  when  solving  a  quadratic  equation  which  detects  surface 
collisions  seems  to  have  eliminated  the  known  occurrences  of  this  effect. 

The  detection  of  multiple  reflections  from  surfaces  and  boundaries  in  an  axisymmetric 
flow  is  a  complicated  task  which  was  addressed  in  a  report  [27]  on  the  2-D  DSMC  code 
which  was  in  use  at  Imperial  College,  London.  The  DSMC2A  code  has  not  been  checked 
against  all  the  possible  errors  mentioned  in  [27]  and  there  may  still  be  some  remaining 
errors  (or  some  new  ones  may  have  been  introduced).  However,  after  extensive  testing  in  an 
enclosed  flow  with  no  weighting  factors,  it  has  become  reasonably  clear  that  any  remaining 
ways  for  particles  to  be  erroneously  deleted  are  fewer  than  before.  The  original  frequency  of 
occurrence  was  small  and  is  smaller  now. 

An  apparent  typographical  error  in  the  main  program  (BFND  for  FND)  where  a  previ¬ 
ously  computed  flow  state  was  read  from  a  file,  has  been  corrected.  The  extent  to  which  the 
errors  mentioned  here  may  be  repeated  in  other  DSMC  codes  supplied  with  [2]  has  not  been 
investigated. 
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At(2RT,)’'*/A, 


Figure  1:  CPU  time  (SPARCstation  10)  required  to  move  119,124  simulator  particles  156 
times  (  tmax  -  156At  ),  sample  the  flowfield  52  times,  and  write  the  final  flow  state  to  disk. 
Gas  at  rest  between  two  coaxial  stationary  cylinders.  For  DSMC  the  non-dimensional  time 
step  shown  on  the  x-axis  is  approximately  equal  to  the  number  of  collisions  per  particle  per 
time  step. 
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diffusely  reflecting  surface  =  T, 


Figure  2:  Computational  domain  and  time  averaged  simulated  flow  showing  toroidal  vortices. 
EPSM.  200x50  grid.  Ntot  =  119, 124.  f  =  41  —  60  revolutions. 
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symmetry  b.c. 


Figure  3;  Kn  =  0.02.  200x50  grid.  EPSM.  Ntot  =  119, 124.  t  =  5  revolutions. 


Figure  4:  Kn  =  0.02.  200x50  grid.  EPSM.  Ntot  =  119,124.  £=10  revolutions. 


Figure  5:  Kn  =  0.02.  200x50  grid.  EPSM.  Not  =  119, 124.  t  =  20  revolutions. 
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Figure  9:  Repeat  run  with  different  random  numbers.  Kn  =  0.02.  200x50  grid.  EPSM. 
Ntot  =  119, 124.  i  =  20  revolutions. 
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Figure  10:  Repeat  run  with  different  random  numbers.  Kn  =  0.02.  200x50  grid.  DSMC. 
Ntot  =  119, 124.  i  =  20  revolutions. 
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Figure  11:  Average  of  23  simulator  particles/cell.  Kn  =  0.02.  200x50  grid.  EPSM.  N-, 
238,248.  f  =  10  revolutions. 


Figure  12:  Average  of  23  simulator  particles/cell.  Kn  =  0.02.  200x50  grid.  EPSM.  Nt 
238,248.  t  =  10  revolutions.  Repeat  run  with  different  random  numbers. 


Figure  13:  Average  of  23  simulator  particles/cell.  Kn  =  0.02.  200x50  grid.  DSMC.  Ntot 
238,248.  <  =  10  revolutions. 


Figure  14:  Average  of  23  simulator  particles/cell.  Kn  =  0.02.  200x50  grid.  DSMC.  Ntot 
238,248.  f  =  10  revolutions.  Repeat  run  with  different  random  numbers. 
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Figure  16:  Kn  =  0.0133.  300x75  grid.  EPSM.  Ntot  =  238,248.  i  =  10  revolutions. 
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Figure  17:  Kn  =  0.0133.  300x75  grid.  EPSM.  Ntot  =  238,248.  i  =  20  revolutions. 
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Figure  18:  Kn  =  0.0133.  300x75  grid.  DSMC.  Ntot  =  238,248.  i  =  5  revolutions. 
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Figure  19:  Kn  =  0.0133.  300x75  grid.  DSMC.  Ntot  =  238,248.  i  =  10  revolutions. 
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Figure  20:  Kn  =  0.0133.  300x75  grid.  DSMC.  Ntot  =  238,248.  £  =  20  revolutions. 
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Molecular  velocity  normal  to  cell  interface  vy{2RT*)^'^ 


Figure  21:  Distribution  of  c„  for  particles  crossing  the  interface  between  two  cells.  Solid 
lines;  distribution  in  EPSM;  dashed  line:  hypothetical  equilibrium  distribution  derived  from 
the  EPSM  distribution 
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Figure  22:  Two  cells  exchange  momentum  giving  rise  to  an  effective  shear  stress  between 
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